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We consider the job assignment problem in a multi-server system 
consisting of N parallel processor sharing servers, categorized into M 
(<S N) different types according to their processing capacity or speed. 

Jobs of random sizes arrive at the system according to a Poisson 
process with rate NX. Upon each arrival, a small number of servers 
from each type is sampled uniformly at random. The job is then 
assigned to one of the sampled servers based on a selection rule. We 
propose two schemes, each corresponding to a specific selection rule 
that aims at reducing the mean sojourn time of jobs in the system. 

We first show that both methods achieve the maximal stability 
region. We then analyze the system operating under the proposed 
schemes as N ^ oo which corresponds to the mean field. Our results 
show that asymptotic independence among servers holds even when 
M is finite and exchangeability holds only within servers of the same 
type. We further establish the existence and uniqueness of stationary 
solution of the mean field and show that the tail distribution of server 
occupancy decays doubly exponentially for each server type. When 
the estimates of arrival rates are not available, the proposed schemes 
offer simpler alternatives to achieving lower mean sojourn time of 
jobs, as shown by our numerical studies. 


1. Introduction. Consider a stream of jobs arriving at a multi-server 
system consisting of a large number of parallel processor sharing servers. 
The servers are categorized into different types or clusters according to their 
processing capabilities. Each job, upon arrival, is assigned to a server where 
it completes its service and leaves the system. The objective is to design job 
assignment schemes that reduce the average sojourn, or response, time of 
jobs in the system. 

1.1. Motivation. The problem of job assignment is central in multi-server 
resource sharing systems that process delay sensitive web requests. Exam¬ 
ples include data centers and web server farms running applications such as 
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online search, social networking etc. In these systems, a small increase in 
the average response time of requests may cause significant loss of revenue 
and users [15]. Therefore, it is critical to reduce the average response time 
of jobs in such systems. 

Reduction in the average response time can be achieved by assigning ar¬ 
rivals to less congested servers [19, 7, 21] in the system. However, in today’s 
systems, where the number of front end servers is large, obtaining state 
information of all the servers incurs a significant communication overhead. 
For such systems, randomized job assignment schemes, in which each assign¬ 
ment decision is made based on comparing the states of a random subset of 
d (> 2) servers, are promising solutions. For systems with identical servers 
(homogeneous), such randomized schemes have been shown [18, 10, 6] to re¬ 
sult in a significant reduction in mean response time of jobs as compared to 
state independent schemes, in which job assignments are made independent 
of server states. This implies that for large homogeneous systems, a small, 
randomly chosen subset of servers is representative of the distribution of 
load in the overall system. 

In this paper, we consider heterogeneous systems where servers are grouped 
into different types or clusters, often geographically separated, based on their 
capacities. Motivated by the aforementioned intuition arising from the ho¬ 
mogeneous case, we consider randomized job assignment schemes, in which 
a small random subset of servers is sampled from each server type. The 
least loaded servers of each type are then compared based on the instanta¬ 
neous processing rates they offer. The job is then assigned to the server that 
provides the highest processing rate. We consider processor sharing (PS) 
as the service discipline in this paper since it closely approximates round- 
robin discipline with small granularity [14] usually employed in server farms. 
Moreover, processor sharing discipline has the desirable property of being 
insensitive to job length distribution type [8]. 

1.2. Related literature. Randomized job assignment schemes have been 
primarily studied in the literature for a system consisting of N identical first 
come first serve (FCFS) servers, which is also referred to as the supermarket 
model. Most studies consider the so called shortest-queue-d (SQ(d)) scheme 
in which each job is assigned to the shortest of d randomly chosen queues. 

For d > 2, [18] showed, using the theory of operator semigroups, that 
the equilibrium queue sizes decay doubly exponentially in the limit as the 
system size increases (as N —>■ oo). Mitzenmacher in [10, 11] derived the 
same result using an extension of Kurtz’s theorem [5]. In [17], a coupling 
argument was used to show that larger values of d result in more even distri- 
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bution of loads among the servers. Chaoticity on path space (or asymptotic 
independence among queue length processes) was established in [6] using 
empirical measures on the path space. Results of [18] were generalized to 
the case of open Jackson networks in [9]. 

Recently, in [3], the SQ(d) scheme was analyzed under more general ser¬ 
vice disciplines and service time distributions. It was shown that in the case 
of FCFS discipline and power-law service time distribution, the equilibrium 
queue sizes decay doubly exponentially, exponentially, or just polynomially, 
depending on the power-law exponent and the number of choices, d. The 
stability of more general randomized schemes for non-idling service disci¬ 
plines was analyzed in [2], which derived a sufficient condition under which 
such networks are stable. Asymptotic independence of servers in equilibrium 
was proposed in [4] under local service disciplines and general service time 
distributions. However, the result was proved only for FCFS service disci¬ 
pline and service time distributions having decreasing hazard rate (DHR) 
functions. 

The tradeoff between sampling cost of servers and the expected sojourn 
time seen by a customer in the supermarket model was studied under a game 
theoretic framework in [22]. It was shown that for arrival rates within the 
stability region of the network, a symmetric Nash equilibrium for identical 
customers exists in which each customer chooses a fixed number of queues 
to sample. 

Recently, in [13, 12], the SQ{d) scheme was considered for a system of 
parallel processor sharing servers with heterogeneous service rates. It was 
shown that, in the heterogeneous setting, random sampling of d servers from 
the entire system reduces the stability region. However, it can be recovered 
using the SQ{d) scheme over a randomly chosen server type. 

1.3. Main results. In this paper, we propose two new randomized schemes 
for job assignment in the heterogeneous scenario. In both the schemes, upon 
arrival of a job, a small number of servers of each type is randomly sampled. 
The sampled servers are then compared based on their states and the arrival 
is assigned to the best server among the chosen set of servers. The metric 
for choosing the best server distinguishes the two schemes. 

This represents a scenario where a centralized dispatcher first requests 
information from each bank or type of servers and then routes the job to 
the server that is going to give the lowest response time among the sampled 
servers. The number of servers sampled from a given type depends on the 
tradeoff between the sampling cost and the likely sojourn time as in the 
supermarket model in [22]. We do not address the precise tradeoffs in this 
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paper suffice to say that we assume that they could be different at each 
server type. We describe the precise mechanisms below. 

In the first scheme, each arrival is assigned to the sampled server with 
the least number of unfinished jobs. In the second, each arrival is assigned 
to the sampled server offering the maximum processing rate per unfinished 
job. Note that, in the both the schemes, the sampled set contains servers 
of all types. We show that such sampling achieves the maximum possible 
stability region. 

We analyze the performance of the proposed schemes in the limit as the 
system size N ^ oo using the mean field approach. Our analysis shows the 
following. 

• The stationary tail distribution of server occupancies decay doubly 
exponentially in the limiting system. We devise indirect methods to 
show this, since, unlike the homogeneous case, closed form solutions 
of the stationary distribution cannot be obtained in the heterogeneous 
scenario. 

• We establish the existence and uniqueness of the equilibrium point 
of the mean field equations in the space of empirical tail measures 
having finite first moment. Our proof, again, differs from the earlier 
works since closed form solutions cannot be obtained. 

• We show that propagation of chaos holds at each finite time and also 
at the equilibrium. In that, we generalize the earlier results on prop¬ 
agation of chaos to systems where exchangeability holds only among 
servers of the same type. 

We also numerically compare the proposed schemes with existing schemes 
for the heterogeneous case. It is observed that the proposed schemes result 
in lower mean response time of jobs in scenarios where arrival rates cannot 
be estimated. 

1.4. Organization. The rest of the paper is organized as follows. In Sec¬ 
tion 2, we describe the system model, the proposed job assignment schemes 
and our notations. We then analyze the proposed schemes in Sections 3, 4, 
and 5. In Section 6, numerical results are presented that compare the schemes 
and determine the accuracy of the theoretical results derived in the paper. 
Finally, we conclude the paper in Section 7 with a summary and a discussion 
on future work. 

2. Model and notations. We consider a multi-server system consist¬ 
ing of N parallel processor sharing (PS) servers. The capacity, C (bits/sec), 
of a server is defined as the time rate at which it processes a single job 


RANDOMIZED JOB ASSIGNMENT SCHEMES 


5 



Fig 1. System consisting of N parallel processor sharing (PS) servers, categorized into M 
types. There are 'fjN servers of type j, each of which has a capacity or rate Cj. Arrivals 
occur according to a Poisson process with rate NX. For each arrival, the job dispatcher 
samples dj servers of type j and routes the arrival to to one of the sampled servers. 


present in it. If there are q{t) jobs present at a server of capacity C at time 
t, then the instantaneous rate at which each job is processed in the server 
is given by C/q{t). Depending on their capacities, the servers in the system 
are categorized into M (<C N) types. Define J = {1, 2,... , M} to be the 
index set of server types. The capacity of type j servers is denoted by Cj, for 
j € J, and we assume, without loss of generality, that the server capacities 
are ordered in the following way: 


(2.1) Cl < ^2 < ... < Cm- 

Further, for each j € J, we denote the proportion of type j servers in the 
system by 7 ^ (0 < 7 ^- < 1). Clearly, 7i = 1- 

Jobs arrive at the system according to a Poisson process with rate N\. 
Each job is of random length, independent and exponentially distributed 
with a finite mean ^ (bits).^ The inter-arrival times and the job lengths are 
assumed to be independent of each other. Upon arrival, a job is assigned to 
one of the N servers where the job stays till the completion of its service 

^As discussed later, our results do not depend on the type of job length distribution 
due to the insensitivity of the processor sharing discipline. 
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after which it leaves the system. The model is illustrated in Figure 1. We 
consider the following two job assignment schemes. 

2.1. Scheme 1. In this scheme, upon arrival of a job, dj servers of type j 
are sampled uniformly at random from the set of N'jj servers of type j, for 
each j G J. Note that this sampling is done at the cluster of type j servers 
by a local router. 

Let ■ ■ •, denote the vector of occupancies of the dj 

sampled servers of type j. For each type j € J", a sampled server with index 
kj is chosen for further comparison where kj is given by 

(2.2) kj = SiTg min 

l<r<dj 1 J 

In case of ties among sampled servers of type j, the index kj is chosen 
uniformly at random from the tied servers of that type. The occupancy 
information of the server corresponding to kj is sent to the central dispatcher. 

Using this information from each of the clusters j ^ J the arriving job 
is assigned by the dispatcher to the type i sampled server having index ki 
where 


(2.3) i = 

Ties across server types are broken by choosing the server type having the 
highest capacity among the tied servers. Thus, in this scheme, each arrival 
is assigned to the server having the least instantaneous occupancy among 
the subset of randomly selected servers. 

2 .2. Scheme 2. As in Scheme 1, upon arrival of a job, a random subset 
of dj servers of type j is chosen uniformly, for each j ^ J. Then from each 
typ6 d ^ J) ^ server with index kj is chosen according to (2.2) for further 
comparison across different server types. The arriving job is finally assigned 
to the type i sampled server having index ki if 

(2.4) i = arg ^ma3^ } . 

Note that the quantity Cjjq^ ^ denotes the processing rate per unfin¬ 
ished job at the sampled type j server with index kj. Thus, in this scheme. 
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an arrival is assigned to the server that provides the highest processing rate 
per job among the sampled set of servers. Ties are broken in the same way 
as described in Scheme 1. 

It is clear that Scheme 2 differs from Scheme 1 only in the criterion for 
server selection. In Scheme 1, server selection is done based only on the in¬ 
stantaneous occupancies of the sampled servers, whereas in Scheme 2 server 
capacities are also taken into account in the selection criterion. Note that 
in the heterogeneous scenario a server with higher occupancy can still pro¬ 
vide a higher processing rate than a server with lower occupancy. Therefore, 
Scheme 2 provides a finer metric for server selection. 

2.3. Notations. We dehne the following real sequence spaces: 

(2.5) = {{gn}n&+ : ffo = l,5n > 5n+i > 0, iVy^c/n € N Vn € Z+}, 

(2.6) U = {{5n}„6z+ : 50 = 1, fi'n > 5n+i > 0 Vn € Z+}, 

OO 

(2.7) U = {{gn}nez+ : 50 = l,5n > 5n+i > 0 Vn € Z+,^ 5 („ < oo}. 

n=0 

Let Yljej^N^ ^ , and denote the Cartesian products of U, and 

U, respectively, over j € An element u = |un \ j £ J,n ^ ^+| belongs 

to Oi'eT 1 1 if 3 € J, the sequence | Un^ | belongs 

to U, or lA, respectively. For u, w € we dehne the following distance 
metric 


( 2 . 8 ) 


|u — w|| = sup sup 
jGjneZ+ 




n -I- 1 


It can be easily verihed that under the metric dehned in (2.8), the space 
is compact (and hence complete and separable). Further, for any k € Z_|_ 
and i,j € N, we dehne 


(2.9) 

( 2 . 10 ) 


lk\ij 




+ 1 ) 


where \x\ denotes the greatest integer not exceeding x and [x] denotes the 
smallest integer greater than or equal to x. 
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Let be a measure space and / : Lf —>• M be a /in-integrable 

function. We define duality brackets as = / fdnH- We denote the 

weak convergence (convergence in distribution) of a sequence of probability 
measures Pn (random variables X„) to a probability measure P (random 
variable X) hy Pn ^ P {Xn X). 

3. Stability analysis. In this section, we derive the sufficient condition 
for the system to have a finite expected number of jobs at all times under 
Scheme 1 and Scheme 2. In other words, we find the set of arrival rates 
for which the Markov process describing the time evolution of the system 
is positive Harris recurrent or stable. We use the stability condition derived 
in [2] for more general join-the-shortest-queue (JSQ) networks. 

Theorem 3.1. The system under consideration is stable under both 
Scheme 1 and Scheme 2 if 

(3.1) XKnJ^nCj- 

HJ 

Proof. Suppose that the N servers in the system are indexed by the set 
S = {1, 2,... , A^}. For each job, we define a selection set to be the subset 
of Yl/jej dj servers sampled at its arrival. We denote by pA the probability 
that the subset T C 5 is chosen as the selection set for an arrival. Note that 
PA, ^ ^ 5, defines the job assignment scheme used. Under Scheme 1 and 
Scheme 2, the probability pA is non-zero only for subsets A which contain dj 
servers of type j for all j € H and for each such a subset A, the probability 
PA is given by 


(3.2) 


PA = 




Now according to Corollary 1.1 of [2], the system under consideration is 
stable if it is subcritical^ i.e., if it satisfies condition (1.2) of [2]. Note that 
the additional conditions (1-11) and (1.12) of Corollary 1.1 of [2] are au¬ 
tomatically satisfied since interarrival times are exponentially distributed. 
Applying condition (1.2) of [2] to the system under consideration, we obtain 
the sufficient condition for stability of the system to be 


-1 


p = max 
BC5 


P ^ C'(n) 
n&B / 


NX'^ Pa\ < I, 


ACIB 


(3.3) 
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where denotes the capacity of the server with index n in the set S. 
Clearly, for Scheme 1 and Scheme 2, the term within the braces in (3.3) is 
non-zero only when the subset B is composed of at least dj servers of type 
j for all j € J. Let Bj (> dj) denote the number of type j servers in B. 
Using (3.2) and (3.3) we now have 


(3.4) 


NX 


p = max 

BCS\N^j>Bj>djy j^J 




(Bj\ 


n 


i&J \d 




BjCj 


- is increasing with respect 


It is easy to verify that that the function 

to Bj for each j € J. Hence, the expression within the braces in (3.4) is 
maximized when we set Bj = N^j. Hence, we have 


(3.5) 


A^A 1 _ A 


Therefore, from (3.3) and (3.5) we conclude that the system under consid¬ 
eration is stable under Scheme 1 and Scheme 2 if (3.1) holds. □ 


Remark 3.1. An alternative proof of stability via a coupling argument 
is as follows: Consider a modified scheme in which, upon arrival of each 
job, one server is chosen from each type uniformly at random (i.e., dj = 1 
for all j G J). The job is then routed to the sampled server of type j 
with probability for each j ^ J. k. coupling argument, similar to 

the one discussed in the proof of Theorem 3 of [9], shows that the system 
operating under the modified scheme always has higher number of unfinished 
jobs than that operating under Scheme 1 or Scheme 2. It is easy to check 
that the system operating under the modified scheme is stable under (3.1). 
Hence, the system operating under Scheme 1 and Scheme 2 also must be 
stable under (3.1). 

As discussed in [2], for A > system under consideration 

is unstable under any job assignment policy. Thus, from Theorem 3.1 we 
conclude that Scheme 1 and Scheme 2 achieve the maximal stability region. 


4. Mean field analysis. We now analyze the time evolution of the 
number of jobs in the system under Scheme 1 and Scheme 2. Its exact 
characterization is difficult since under both the schemes, arrivals at a given 
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server depend on the states of other servers. However, it is possible to analyze 
the system in the limit as the system size N ^ oo. Such a limit is known as 
the mean field limit [10, 18, 9] and it exists because under random sampling 
of a fixed number of servers from each type the statistical properties of 
the system do not change when states among servers of the same type are 
permuted. 

To formally state our results, we define the process 

(4.1) XAr(t) = j e J,nG Z+l for t > 0, 

where denotes the fraction of type j servers having at least n un¬ 

finished jobs at time t. Thus, |x^^^(t),n € Z+| denotes the empirical tail 
distribution of occupancy of type j servers at time t. Clearly, XAr(t) is a 
Markov process in the state space n, 

4.1. Convergence to the mean field. The main aim of this subsection is 
to prove the following result. 

Theorem 4.1. //x7v(0) converges in distribution to some constant g G 

as N ^ oo, then the process converges in distribution to 

a process {u(t)}j^g, lying in the space as N ^ oo. For Scheme 1, the 
process u(t) is given by the solution of the following system of differential 
equations 


(4.2) u(0) = g, 

(4.3) u{t) = l(u(f)), 

where the mapping 1 : —)■ (M^+)’^ is given by 


(4.4) 

(4.5) 


l^\n) = 0, for j G J, 




X f /' (j) \dj 


u 


k-l 




i-1 


M 


i=l 


w TJ 


n (4A) n 

%=j+i 


- fiCj , for k> l,j G J. 


For Scheme 2, the process u(t) is given by the solution of 
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(4.6) u(0) = g, 

(4.7) u{t) = l{u{t)), 

where the mapping 1 : —>■ (M^+)’^ is given by 


(4.8) 

(4.9) 


S^(u) = 0, for j € J, 


?’(u) = - 

7J 

M 


U 


(i) 

fc-i 


- <nf)‘ 


J-l 

n 

i=l 




di 


di 


^ n - 4+i) , fork >1J ej. 

i=j+l 


The process {u(i)}^>Q, defined in the theorem above, is referred to as 
the mean field. We first note that Theorem 4.1 implicitly assumes that the 
ordinary differential systems (4.2)-(4.3) and (4.6)-(4.7) have unique solutions 
in the space . In the following proposition, we show that this is indeed 
the case. To emphasize the dependence of the solution u(t) on the initial 
point g, we will often denote u(t) by u(t, g). 

Proposition 4.1. // g G , then each of the systems (4.2)-(4.3) 

and (4.6)-(4.7) has a unique solution u(t, g) G , for all t > 0. 


Proof. The proof is given in Appendix A. 


□ 


We will prove Theorem 4.1 using the theory of semigroup operators of 
Markov processes as in [18, 9]. Before doing so, we recall the following 
from [5]. 

• For the process {xAr(t)}^>Q, the operator semigroup {TAr(t)}^>Q acting 

on continuous functions / : —7> M is defined as 

TAr(t)/(x)=E[/(xAr(t))|xAr(0)=x] Vt>0,xG 

j&J 

• For the deterministic process {u(t)}^>g, the transition semigroup {T(t)}j>Q 
acting on continuous functions / : hl^ ^ R is defined as 


T(t)/(x) =/(u(t,x)) Vt>0,xGW^. 
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In the next proposition, we show that TAr(t) converges to T(t) uniformly 
on bounded intervals. This in conjunction with Theorem 2.11 of Chapter 4 
of [5] proves Theorem 4.1. 

Proposition 4.2. For both Scheme 1 and Scheme 2, and for any con¬ 
tinuous function f : M and t > 0, 


(4.10) 


lim 

A—>00 


sup 


TA(t)/(g) - /(u(t,g))| = 0 


and the convergence is uniform in t within any bounded interval. 


Proof. The proof is given in Appendix B. 


□ 


Remark 4.1. We note that Theorem 4.1 implies that if X7v(0) => g € 
as —>■ 00 , then the following weaker convergence results also hold: 

1. For each t > 0, xjv(t) ^ u(t, g) as N ^ 00 . 

2. For each t > 0, j ^ J, and k € x^^\{t) u^^\t, g) as —>■ 00 . 


3. For each t > 0, j G and k € E 


X 


U) 


g) as > 00 . 


The last assertion follows from the first since x''^ i^{t) is bounded for each 
N,j,k,t. 


4.2. Properties of the mean field. In this section, we characterize some 
important properties of the mean field. In particular, we show that, under 
the stability condition (3.1), both (4.2)-(4.3) and (4.6)-(4.7) have unique 
equilibrium points in . Further, we show that the equilibrium points are 
globally asymptotically stable for both systems. 

Let P, P denote the equilibrium points of (4.2)-(4.3) and (4.6)-(4.7), 
respectively. In other words, P and P satisfy 1(P) = 0 and 1(P) = 0. Hence, 
for all k G Z+ and f G J the following must hold 


(4.11) 


p(i) 

^k-ei 


p(i) 

^k+2 


= A, 



di'o'O 




X n (T' 

) ■ n (k 

i=l 

i=j+i 


W 


di 
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(4.12) 


p(i) _ p(i) 

^k+2 


= A, 



■pis) VA 


A+4 j 




xn(A'i. 

in 

i=l 

i=j+i 


P, 


(i) 


A (?) 

where Aj = for each j G J. Note that by definition we have Pq = 

Pg = 1 for all j € J. The next proposition reveals an important property 
of the equilibrium points P and P. To state it we first need the following 
definition. 


Definition 4.1. A real sequence {zn}n>i said to decrease doubly ex¬ 
ponentially if and only if there exist positive constants L, oj < 1, 9 > 1, and 
K such that Zn < kw®" for all n> L. 

Hence, if a sequence {zn]n>i decays doubly exponentially, then it is 
summable, i.e., < oo. 

Proposition 4.3. Assume that for each j G J, J- 0 ^ ^ oo. 

Then the following equations must hold 


(4.13) 


E 

j&J 




3&J 


(4.14) 


P 


( 1 ) 


/+! 


M P\ 

+E' 

i=2 


(j) 

^ fp(i) 


A, 


di 


M 


n T, 


■Ai) 

U-iJi 


i=2 


Further, for each j G J, the sequences ^Plf\k G Z+| and ^P^\k G Z+| 
decrease doubly exponentially. In particular, under the assumption of the 
proposition, both | P^^\ k G and | P^^\ k G Z+| are summable sequences. 


Proof. We prove the proposition for P. The proof for P follows along 
the same line of arguments. For a fix j we add (4.11) for all /c > Z and use 

fj\ 

limfc_^oo = 0 to obtain 






14 


A. MUKHOPADHYAY, A. KARTHIK, AND R. R. MAZUMDAR 


(4.15) 
)(i) 


pli\ = A E 


k>l 


n (p, 


di 


M 


y{i) 

k+1 


di 


i=l 


i=j+l 


i-1 

n 

2=1 


,{*) 


di 


M 


HR, 


fc+i 


di 


l=J 


Now, multiplying both sides of the above equation by and adding over 
all j ^ J and using liuik^ooPjf^ = 0 yields (4.13). From (4.13) we obtain 
< Wjej {^k^) ' ^ (-Pfc) > where Pk = maxi<j<M Pk^ and d = 

Thus, we have Pk^i < SPk, where 6 = (^Pk'^ niaxi<j<M(Aj). 

(i) 

Since by hypothesis, for each j, ^ —)> 0 as /c —>■ oo, one can choose k 
sufficiently large such that J < 1. Hence, we have Pk\'^ < 

6Pk. Similarly we have, (iiiaxi<j<Af ) < d Pk- This proves that the 
sequence | Pk\ k € Z+| decreases doubly exponentially for each j. □ 


The following proposition guarantees that there exists equilibrium points 
of systems (4.2)-(4.3) and (4.6)-(4.7) in . 


Theorem 4.2. Under condition (3.1), there exists an equilibrium point 
P of the system (4.2)-(4.3) andP of the system (4.6)-(4.7) in the space . 


Proof. The proof is given in Appendix C. □ 

The next theorem shows that P and P are the unique globally asymptot¬ 
ically stable equilibrium points of the systems (4.2)-(4.3) and (4.6)-(4.7) in 
the space lA^. 


Theorem 4.3. Under condition (3.1), 


(4.16) lim u(t, g) = P € for all g G , 

t^OO 

for Scheme 1 and 


(4.17) lim u{t, g) = P G for all g G 

t^OO 

for Scheme 2. Hence, P and P are globally asymptotically stable fixed points 
of systems (4.2)-(4.3) and (4.6)-(4.7), respectively. Furthermore, P and P 
are the only equilibrium points of the above systems in the space . 
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Proof. The proof for Scheme 1 is given in Appendix D. For Scheme 2, 
the theorem can be similarly proved. □ 

We now show that, under (3.1), the stationary distribution of the process 
XAT converges weakly to the Dirac measure concentrated at the unique equi¬ 
librium point of the mean field. Let tin denote the stationary distribution 
of the process x^v- Clearly, ttat exists and is unique under (3.1). Further, 
for each fixed N, X 7 v(t) ^ XAr(oo) as t —>■ oo, where xn{oo) is a random 
variable distributed as ttat. 

Theorem 4.4. Under condition (3.1), we have 


(4.18) 

for Scheme 1 and 


(4.19) ttn ^ dp, 

for Scheme 2. 


Proof. We prove the theorem for Scheme 1. The proof for Scheme 2 
follows similarly. 

Note that since the space is compact, so is the space of probability 
measures on . Therefore, the sequence of probability measures 
has limit points. Thus, in order to prove the theorem we need to show that 
all limit points coincide with Jp. 

Due to Theorem 4.1, any limit point tt of the sequence tt^v must be 
an invariant distribution of the maps g i-A u(t,g). Hence, by uniqueness 
proved in Theorem 4.3, it is sufficient to prove that tt is concentrated on 
lA^. To prove that vr is concentrated on lA^ it is sufficient to show that 




Z^n>l yn 


implies that En>i5n ^ 


E. 




(i) 


TTiV 9n 


< OO for all j G The coupling described in Remark 3.1 

^ <1. Hence, 

□ 


< where p = 




E, 




< This completes the proof. 


We have so far established that the interchange property indicated in 
Figure 2 holds. Note that the convergences indicated in the figure are in 
distribution. 
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Xjv(t) 


Xjv(oo) 


N ■ 


oo 


Theorem 4.1 


Theorem 4.4 
N ^ oo 


u{t) 


Fig 2. Commutativity of limits 


4.3. Propagation of chaos. In this subsection, we focus on the occupan¬ 
cies of a given finite set of servers as N —> oo. We show that as the system 
size grows the server occupancies become independent of each other. Such 
independence holds at any finite time and also at the equilibrium, provided 
that the initial server occupancies satisfy certain assumptions. This is for¬ 
mally known as the propagation of chaos [6, 16] or asymptotic independence 
property [4, 3] in the literature. 

To formally state the results we introduce the following notations. Let 
q^]^’^\t), for 3 C J and k G {1,2,... ,A^ 7 j}, denote the occupancy of the 
^th ggpygj. Qf type 3 at time t > 0. By q^^’^\oo) we denote the occupancy of 

the server of type j in equilibrium. Further, let 3 ^ ff and 

n G Z_|_, denote the fraction of type j servers having occupancy n at time t > 

0. Dehne the process XN{t) = {x^N,ni^)^3 ^ J,nC Z+|. Clearly, x%\t) = 
G Z+l denotes the empirical distribution of occupancies of type 
j servers and for each n,j, we have XN]n(^) = ® 7 v!n(^) “ 

X%\oo) we will denote the empirical distribution occupancies for type j 
servers in equilibrium. Let the process Q(t) = ^Qn\t),j & J,n G Z+l be 
defined as Qn\t) = Un\t) — ul^]_i{t), for t G [0,oo]. Further, we denote by 
Q^^\t) the distribution on Z+ given by Q^^\t) = |(3n\n G Z+|. We also 
define the following notion of exchangeable random variables. 

Definition 4.2. Let 1 < ^ ^ —3 ^ denote a collec¬ 

tion of N random variables among which Nfj belong to a particular class j 
and are indexed by k, where 1 < k < N'jj. The collection is called intra-class 
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exchangeable if the joint law of the collection is invariant under permutation 
of indices, 1 <k< N'jj, of random variables belonging to the same class. 

Proposition 4.4. For the model considered in this paper, for both schemes, 
1 < ^ 1 < j < -^1 intra-class exchangeable and i/xAr(0) 

g G as N ^ oo, then the following holds 

(i) For each fix k and t G [0,oo], q^j^’^\t) U^^\t) as N ^ oo, where 

U^^\t) is a random variable with distribution Q^^^t). 

(a) Fix positive integers ri, r 2 ,..., vm- For each t G [0, oo], 

l<k<rj,l<j<M^^ I < k < rj,! < j < , 

as N ^ oo, where 1 < k < rj,l < j < M, are independent 

random variables with having distribution for all 1 < 

k <rj. 

Proof. Note that the first part of the proposition is a special case of the 
second part. Hence, it is sufficient to prove the second part. We will provide 
a proof for the M = 2 case. The proof can be readily generalized to any 
M > 2. 

Due to the dynamics of the system (under Scheme 1 or Scheme 2) and 
the hypothesis of the proposition {q^^’^\t),l < k < < j < M} is 

intra-class exchangeable for all t G [0, oo]. The hypothesis of the proposition 
also implies that XNi't) ^ Q(t) as N" —)• oo for all t G [0, oo]. Henceforth, we 
will omit the variable t in our calculations, which hold for all t G [0, oo]. 

To prove the proposition, it is sufficient to show that the following con¬ 
vergence holds as —)■ oo. 


(4.20) 


E 


ri r2 


ri r2 


k=l k=l 


for all bounded mappings ^ 1^+- Now we have 
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(4.21) 


E 


ri 


n wn n 


r2 


-xm 


.k=l 


k=l 


ri 


r2 


< 


E 


ri 




r2 


xm 


Lfc=l 


k=l 


-E 


+ 


E 


ri 


) Yli^k,X^N) 


f2 


( 2 )\ 


.k=l 


k=l 


k=l k=l 

fi{(pk,xi'’) n(V'fc,xw) 

.k=l k=l 

ri r2 


k=l 


k=l 


Note that the second term on the right hand side of the above inequality 
vanishes as —>■ oo since x^n ^ as N ^ oo for j = 1,2 and and 
are constants. Now, due to exchangeability we have 


(4.22) E 




42A) 


Lfc=l 


k=l 


(iV7i)n(iV72) 


r2 


X E 


ri 


E E 

uGP{ri,N'yi) cr'GP{ri,N''yi) k=l k=l 


r2 


{2,a'{k)) 


where {N)^ = N{N — 1)... {N — k + 1), and P{r,n) denotes the set of 
all permutations of the numbers {1, 2,... , A^} taken r at a time. Also, by 
definition of Xn have 


(4.23) 


E 


ri 


J’2 




= E 



1 

A^7i 


A 71 

1=1 


( 1.0 

Qn 


X 



Hence, the first term on the right hand side of (4.21) can be bounded as 
follows 


E 


n^^(*‘’)n^^(«r) 


- E 


ri 


r2 


_k=l k=l 


< ( 1 - 


(lV7i),,(iV72) 


r2 


{NjiY^iNx^y^ 


0 as Ai —>■ 00 , 
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where max (||(/)fc||oo, HV’fcHoo) = B. This completes the proof. □ 

Thus, the above proposition shows that in the limiting system server 
occupancies become independent of each other. It also shows that the sta¬ 
tionary occupancy distribution of any type j server is given by Q^^\oo) = 
G Z+l for Scheme 1 and Q^^\oo) = ^Pn'^ — P^}_i,n G Z+| 

for Scheme 2. 

5. Computation of the stationary distribution. So far we have 
shown that in the limiting system {N oo) each finite collection of servers 
behave independently and the stationary tail distribution of occupancy of a 
typ6 J G server in the limiting system is given by ^P^\k G Z+| under 

Scheme 1 and G under Scheme 2. Using the independence of 

servers in the limiting system we conclude the following proposition. 

Proposition 5.1. In equilibrium, the arrival process of jobs at any given 
server in the limiting system is a state dependent Poisson process. Further, 
the arrival rate of jobs to a server of type j € H when it has occupancy k in 
the equilibrium is given by 


(5.1) 


= 


P, 


U) 


- P 


,(i) 


fc+i 


7j 


Pu^ - P, 


U) 

k+l 


n 

2=1 




M 


n 

i=j+i 


p 


(i) 


k+l 


di 


for Scheme 1 and 


(5.2) 


= 


P, 


U) 


- P 


yU) 


k+l 


7i 


^ - P, 


U) 


k+l 


i-1 

n 

2=1 


Pr 


(7 

W,' 


M 


n 

i=j+i 


p 


(d 


for Scheme 2. 


Proof. We provide the proof for Scheme 1. The proof for Scheme 2 
follows from similar line of arguments. 

Consider a tagged type j server in the system and the arrivals that have 
the tagged server as one of its possible destinations. These arrivals constitute 
the potential arrival process at the tagged server. The probability that the 
tagged server is selected as a potential destination server for a new arrival is 
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= 7 v^- Thus, due to Poisson thinning, the potential arrival process 
( dp ^ 

to the tagged server is a Poisson process with rate x NX = 

l\'yj 7j 


Next, we consider the arrivals that actually join the tagged server. These 
arrivals constitute the actual arrival process at the server. For finite N, 
this process is not Poisson since a potential arrival to the tagged server 
actually joins the server depending on the number of jobs present at the other 
possible destination servers. However, as N ^ oo, due to the asymptotic 
independence property shown in 4.4 the occupancies of the sampled servers 
become independent of each other. As a result, in equilibrium the actual 
arrival process converges to a state dependent Poisson process as N ^ oo. 

Consider the potential arrivals that occur to the tagged server when its 
occupancy is k. This arrival actually joins the tagged server with probability 
when X other servers among the dj servers of type j have occupancy 
k, all the di servers of type i < j have at least occupancy k, and all the d* 
servers of type i > j have at least occupancy k + 1. Thus, the total arrival 
rate Xp can be computed as 


(5,3) Af 


7i 



dj 


X 


{N - hh)' 

(hh)" 

— l—x 




X n (N 

)■ n 

( 

7 fc+i 

i=l 

z=j+l 



which simplifies to (5.1). □ 

Hence, the above proposition shows that in equilibrium the arrival rate 
at a given server depends on the stationary tail probabilities PP, k € Z_|_ 
and j € J. 

The stationary tail probabilities can in turn be expressed as functions of 
the arrival rate. Indeed, in equilibrium the global balance equations (which 
hold under state dependent Poisson arrivals due to Theorems 3.10 and 3.14 
of [8]) yield 


(5.4) for j £ J,k 

where — Pjfji- Hence, the equilibrium point P is the unique fixed 

point of the mapping 0 : defined as 0(P) = F{G(P)), where G(-) 
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denotes the mapping from lA^ to the space of possible arrival rates (defined 
by (5.1)) and F{-) denotes the mapping from the space of possible arrival 
rates to the space (defined by (5.4)). Thus, the equilibrium point P can 
be computed using the fixed point iterations (i.e., by repeatedly applying 
the mapping ©(•) to some arbitrary point Q G .) 

Remark 5.1. So far our results have been obtained for exponential job 
length distributions. Note that the conclusions of Proposition 5.1 continue to 
hold for any job length distributions due to the Whittle balance criterion [20] 
that can be shown to hold for the stationary distribution (also see Theorems 
3.10 and 3.14 of [8]). In view of the uniqueness of the stationary distribution 
and propagation of chaos this suggests that in stationarity the servers are 
asymptotically independent for general job size distributions. In Section 6, 
we provide numerical evidence to support insensitivity. 

Remark 5.2. From Proposition 4.4 it directly follows that the expected 
occupancy of a type j server at equilibrium is given by Scheme 1 

and Scheme 2. Hence, a simple application of the Little’s law, 

yields that the mean sojourn time of jobs in the limiting system is given 

by f = 1E" 1 Er.i liPt’ for Scheme 1 and f = ^ Efli E^.l for 

Scheme 2. Thus, the mean sojourn time of jobs in the limiting system can be 
computed using stationary tail probabilities which in turn can be computed 
using the fixed point method described in this section. 

6. Numerical Results. In this section, we present simulation results 
to compare the different job assignment schemes discussed in this paper. The 
results also indicate the accuracy of the asymptotic analyses of the Scheme 1 
and Scheme 2 in predicting their performance in a finite system of servers. 
We set /r = 1 in all our simulations. 

To determine accuracy of the asymptotic analysis presented in the paper 
we first compare the results obtained from the theoretical analysis with that 
obtained from the simulations. In Figure 3, we plot the mean sojourn time 
jobs as a function of the normalized arrival rate. A, for different values of 
the system size N. We observe a very good match between the analysis 
and simulation results for N = 100. For N = 10 and N = 20 the relative 
errors between the analysis and the simulation results are around 10% and 
5%, respectively. Thus, we conclude that the asymptotic analysis accurately 
captures the behaviour of the system for moderately large system sizes. 

We now compare the performance of the proposed schemes with that of 
other existing schemes for heterogeneous scenario. In particular, we consider 
the following two schemes as benchmarks. 
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Fig 3. Mean sojourn time jobs as a funetion of X for different values of N. We set 
Cl = 2/3, C 2 = 4/3, and 71 = 72 = 0.5. 



Fig 4. Mean sojourn time jobs as a function of A for different schemes. We set M = 2, 
Cl = 1/5, C 2 ~ 9/5, 71 = 72 = 0.5, and di = d 2 = 2 . Routing probabilities for the state 
independent scheme and the Hybrid SQ(d) scheme are optimized based on A. 
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6.1. The state independent scheme. As a baseline, we consider a scheme 
that assigns an incoming job to a server with a fixed probability, independent 
of the current state of the servers in the system [1]. We denote by pj, for 
j G J, the probability with which an arrival is assigned to one of the servers 
of type j. The probabilities pj, j G J, can be chosen chosen such that 
the mean sojourn time of the jobs is minimized. Clearly, in this scheme, no 
communication is required between the job dispatcher and the servers as the 
job assignment decisions are made independently of the state of the servers. 

6.2. The hybrid SQ(d) scheme. In this scheme [13], upon arrival of a 
new job, the router first chooses a server type j & J with probability pj. 
Then dj servers of type j are chosen uniformly at random from set of N'jj 
servers of type j. The job is then assigned to the server having the least 
number of unfinished jobs among the dj chosen servers. Ties are broken by 
tossing a fair coin. As in the state independent scheme, the probabilities pj, 
j G J, can be chosen such that the mean sojourn time of jobs in the system 
is minimized. 

We choose the parameter values as follows: M = 2, Ci = 1/5, C 2 = 9/5, 
7 i = 72 = 0.5, and di = d 2 = 2. Under this parameter setting, the stability 
region for all the schemes under consideration is A < 1. In Figure 4, we plot 
the mean sojourn time of jobs as a function of the normalized arrival rate, 
A, for Scheme 1, Scheme 2, the state independent scheme, and the hybrid 
SQ(d) scheme. We choose the optimal routing probabilities pj, j G ^7, for 
both state independent scheme and the hybrid SQ,{d) scheme. We observe 
that the mean sojourn time of jobs under Scheme 1 and is almost the same 
as that under Scheme 2 for small values of A. However, for larger values of A, 
Scheme 2 outperforms Scheme 1. This is expected for reasons explained in 
Section 2. We also see that hybrid SQ(d) scheme results in a smaller mean 
sojourn time of jobs than that in Scheme 1 and Scheme 2, for smaller values 
of A. This is because, in the hybrid SQ(4) scheme, the routing probabilities 
are chosen optimally based on the arrival rate A. However, for larger values 
of A, we observe that Scheme 2 outperforms the hybrid SQ(d) scheme. 

To observe the effect of fixing the routing probabilities for the hybrid 
SQ(d) scheme and the state independent scheme, we choose pi = ^ 

for each server type i ^ J. This choice of routing probabilities ensures that 
all arrival rates in the maximal stability region can be supported by the 
system operating under either the state independent scheme or the Hybrid 
SQ(d) scheme. We choose the same parameter setting as before and plot 
mean sojourn time of jobs as a function of A in Figure 5 for the schemes un¬ 
der consideration. In this case, we notice that both Scheme 1 and Scheme 2 
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Fig 5. Mean sojourn time jobs as a function of X for different values of N. We set M = 2, 
Cl = 1/5, C 2 = 9/5, 71 = 72 = 0.5, and di — d 2 = 2. Routing probabilities for the state 
independent scheme and the hybrod SQ(d) scheme are not optimized. 


outperform the hybrid SQ,{d) scheme. Hence, in the scenarios where esti¬ 
mation of arrival rates is not possible, Scheme 2 is a better choice than the 
hybrid SQ,{d) scheme. 


Table 1 

Insensitivity of Scheme 1 


A 

Mean sojourn time T 

Constant 

Power Law 

(Theoretical) 

(Simulation) 

(Simulation) 

0.2 

0.8076 

0.8106 

0.8098 

0.3 

0.8609 

0.8642 

0.8640 

0.5 

0.9809 

0.9852 

0.9840 

0.7 

1.1696 

1.1759 

1.1757 

0.8 

1.3687 

1.3741 

1.3740 

0.9 

1.7531 

1.7641 

1.7645 


We now numerically investigate the behaviour of the proposed schemes 
under different job length distributions. In Table 1, mean sojourn time of jobs 
under Scheme 1 is shown as a function of A, for the following distributions. 

1. Constant: We consider job length distribution having the cumulative 
distribution given by F{x) = 0 for 0 < x < 1, and F{x) = 1, otherwise. 
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2. Power law: We consider job length distribution having cumulative dis¬ 
tribution function given by F{x) = 1 — l/4x^ for x > ^ and F{x) = 0, 
otherwise. 

For both distributions we have = 1. We choose the following parameter 
values M = 2, Cl = 4/3, C 2 = 2/3, N = 100, 71 = 72 = 5 , and di = d 2 = 2. 
We observe that there is insignificant change in the mean sojourn time of 
jobs when the job length distribution type is changed. The results, therefore, 
justify the insensitivity property as discussed in Remark 5.1. 

7. Conclusion. We considered randomized job assignment schemes in 
a multi-server system consisting of N parallel processor sharing servers, cate¬ 
gorized into M (-C N) different types according to their processing capacity 
or speed. In the proposed schemes, a small number of servers from each 
type is sampled uniformly at random at each arrival instant. It was shown 
that due to such sampling the schemes achieve the maximal stability region. 
Mean field analysis was carried out to show that asymptotic independence 
among servers holds even when M is finite and exchangeability holds only 
within servers of the same type. The existence and uniqueness of stationary 
solution of the mean field and doubly exponentially decreasing nature of the 
tail distribution of the number of jobs was established. Numerical studies 
have shown that, when the estimates of arrival rates are not available, the 
proposed schemes offer simpler alternatives to achieving lower mean sojourn 
time of jobs. 


APPENDIX A 

We will prove Proposition 4.1 only for the system (4.2)-(4.3). The proof 
for the system (4.6)-(4.7) follows similarly. 

Define 6{x) = [min(x, I)]+, where [z]^ = max{ 0 , 2 ;} and let us consider 
the following modification of (4.2)-(4.3): 

(A.I) u(0) = g, 

(A.2) u{t) = i(u(f)), 

where the mapping 1 : —>■ ( 1 ^^+)^^ is given by 
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(A.3) i«(u) 
(A,4) ijfhu) 

M 

X n 

i=j+i 


= 0, for j G J, 


\ 

7j L 


e{u 


,(i) 

‘‘k-i 


-{e{u 


,U) 




i=l 


e u 




^Cj 


0 ( u\p \-d{u 


,(i) 

‘‘k+l 


, for A: > 1, j G J. 


Clearly, the right hand side of ( 4 . 5 ) and (A. 4 ) are equal if u G . Therefore, 
the two systems must have identical solutions in . Also if g G , then 
any solution of the modified system remains within . This is because 
of the facts that if Un\t) = for some j, n, t, then /n\u(t)) > 0 

and ip^^{u{t)) < 0, and if Un\t) = 0 for some j, n, t, then lP\u{t)) > 0. 

Hence, to prove the uniqueness of solution of ( 4 . 2 )-( 4 . 3 ), we need to show 
that the modified system (A. 1 )-(A. 2 ) has a unique solution in (]R^+)^. We 
now extend the distance metric defined in (2.8) to the space (M^+)^. 

Using the metric defined in ( 2 . 8 ) and the facts that \x+ — y+\ < \x — y \ for 
any x, y G M, \aib^ — 026^] < |ai — 02! + m |6i — 62I for any ai, 02, 61,62 G 

[ 0 , 1 ], and \ 9 {x) — 9 {y)\ < |x — y| for any x,y G K we obtain 


(A. 5 ) l|l(u)||<i^i, 

(A.6 ) I|i(u) - i(w)|| < A' 2 ||u - w||, 

where u, w G (M^+)^, Ki and K2 are constants defined as Ki = + 

/i(maxj£j-Cj) and K2 = + 3 //(maxi<j<M Cj). The uniqueness 

now follows from inequalities (A. 5 ) and (A.6) by using Picard’s iteration 
technique since (M^+)^ is complete under the metric defined in (2.8). □ 

APPENDIX B 

We prove Proposition 4.2 by showing that the generators of the corre¬ 
sponding semigroups converge as N —>■ 00. We first recollect the following 
from [ 5 ]. 

• The generator of the semigroup {Tjv(t)}i>o acting on functions 
/ : ^ M is given by A^fig) = Eh^g 9gh (/(h) - /(g)), 

where ^gh, with g,h G ^ denotes the transition rate from 

state g to state h. 
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• The generator A of the semigroup {T(t)}j>Q acting on functions / : 

R having bounded partial derivatives is given by A/(g) = 

In the following lemma, we characterize the the generator Ajv associated 
with the process XAr(t). 

Lemma B.l. Let g € be any state of the process XAr(t) and 

e(n, j) = (el!M be the unit vector with Cn^ = 1 and =0ifi^j 

\ ^ J kGZ+,iej ^ 

and k n. Under Scheme 1, the generator Ajv of the Markov process X]\f{t) 
aeting on funetions f : Ojli ^ R is given by 


M 


(B.l) ANf{g) = NX^Yl 

j=l n>l 


M) 

hn-l 


- 


n i-1 


2=1 


n 


di 


M 


fn 


X n [ 9 ^ 

i=j+l 
M 

+ pN EE Cj 

n>lj=l 


di 


/(g+ )-/(g) 




qU) _ o(j) 
iJn l/n+l 


ry e(n,j) 

/(g-^^)-/(g) 


Under Scheme 2, the generator of the Markov process x^(t) acting on 
functions f : ^ R is given by 


M 


(B.2) ANf{g) = NXj2Yl 

j=l n>l 


M) 

i/n-l 




j-1 


,(d 


n (^rn-iL, 


di 


i=l 


M 






i=j+l 


/(g + 


Ay, ■ 


-/(g) 


M 


n>lJ=1 


ffn ^ - ^i+l 


e(n,j) 

/ g - ^TT- - /(g) 


Proof. We only prove the lemma for Scheme 1. For Scheme 2, it can be 
shown similarly. 

We first consider an arrival joining a server of type j with exactly n — 1 un¬ 
finished jobs, when the state of the system is g. This corresponds to the tran¬ 
sition from state g to the state g -|- The term 
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" 1 ^ ]\/^ f ^^ 

^U.iZl[9n-l) Hi =j+i \9n ) denotes the probability with which an ar¬ 
rival joins a type j server with exactly n — 1 jobs. This is because a job joins 
a server of type j with exactly n — 1 occupancy only when the following 
conditions are satisfied: 

• Among the dj sampled servers of type j, at least one has exactly n — 1 
jobs and the rest of them have at least n jobs. 

• For each i < j, all the di sampled servers of type i have at least n — 1 
jobs. 

• For each i > j, all the di servers of type i have at least n jobs. 

Since the arrival rate of jobs is NX, the rate of the above transition is given 
by 


(B.3) 


"g.g+ 


N-tj 


NX 


eN - bp'"' 



2=1 i=j-\-l 


Further, the rate at which jobs depart from a server of type j having exactly 
n jobs is fiCjNjj ■ The expression (B.l) now follows directly 

from the definition of Atv. Q 


We now show that the solutions u(t,g) of (4.2)-(4.3) and (4.6)-(4.7) are 
smooth with respect to the initial point g and their partial derivatives are 
bounded. 


Lemma B.2. For each j, n, j', n', i, k, and t > 0, the partial derivatives 

exist for g € and satisfy 


(B.4) 

and 

(B.5) 




dg. 


(j) 


924*^ (Lg) 


d‘^u{’{t,g) 

dg^^" 

•) 



< exp(Bit) 


< ^(exp(2fiit) - exp(Bit)), 


^dere Cj), and B 2 = 


Proof. The proof follows the same line of arguments as the proof of 
Lemma 3.2 of [9]. We omit the details. □ 
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Proof of Proposition 4.2. The proof is essentially the same as the 
proof Theorem 2 of [9]. We omit the details. □ 

APPENDIX C 

We prove the existence of equilibrium point for Scheme 1. Similar argu¬ 
ments apply for Scheme 2. For simplicity of exposition, we further restrict 
ourselves to the M = 2 case. However, the proof can be extended to any 
M >2. 

The idea is to construct sequences k € Z_|_| for j = 1,2 such that 

they satisfy the following three properties 

P.l Equation (4.11) for j = 1,2. 

P.2 > Pi^\ > 0 for all k G Z+, j = 1,2. 


P.3 P, 


(j) 


0 as A: —7> oo for j = 1,2. 


According to Proposition 4.3, we see that P = ^pj^^\k G Z_|_,j G {1,2}| 

with components Pf^ ' satisfying the above properties, must be an equilib¬ 
rium point of the system (4.2)-(4.3) and also must lie in the space Note 
that if (P.l) holds and pjf'^ > 0 for all k and j, then pjf^ > 

We now construct the sequences |p/^^(a),Z G and |p/^\q;),Z G 
as functions of the real variable a as follows: 


(C.l) 

Pi^\a) 

= 1. 

(C.2) 

pI>"\») 

= 1. 

(C.3) 

pW(«) 

= a. 

(C.4) 

P?\») 

= Aa ^ 

(C.5) 


_ p{i) 

“ P+i 

(C.6) 


_ p{2) 

“ P+l 


a 




(a) - Ai 

X (^S(«)) 

(a) - Aa ^(^P/^^(a)) 


- P, 


,{i) 


i+i 


(«)) 


di 


d2 


,z > 0 


di 


X (p«(a)) \l>0 


Combining the above relations we obtain 
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(C.7) 




n , for Z > 0 

i=i 


Note that that the sequences G Z+| and G Z_|_| are 

constructed such that they satisfy property (P-1). Hence, the the proof will 
be complete if for some a G (0,1) the properties (P.2) and (P.3) are satisfied. 
We first proceed to find a G (0,1) such that the sequences |p/^^(a), / G Z+| 

and |p/^^(q:), I G Z+| are both positive sequences of real numbers in [0,1]. 
This will ensure that (P.2) is satisfied. 

Note that p/^^(l) = 1 for all I G Z+. Hence, from (C.4) we have P^^^^(l) = 
A 2 (1 — and from (C.6) we have 


(C.8) Pfllil) = P® (1) - A 2 ((P/')(1))"' - (pS( 1))"') for / > 0 
Notice that the stability condition (3.1) reduces to 


(C.9) 


1 1 


which implies that P^ Tl) < 1. We claim that there exists some I > 1 such 

{'2'] (' 2 ') 

that P; pi) < 0. Let us assume this is not true. Therefore, P^^ pi) > 0 for 

all I > 0. By (C.8), this implies that |p/^^(l),Z > o| is a non-decreasing 

sequence of numbers in [0,1). Hence by monotone convergence theorem 
(2) 

lim;_>.oo P; rI) exists. Let this limit be denoted by /3, where 0 < /3 < 1. 

f 2^ 

Thus, adding (C.8) for / > 0 and using lim;_>.oo P; pi) = /3 we obtain 


1 - 


Ai 


1^2 


>/3 


1 -^] + 1 -/ 3 '"^ 


/ I \ ^_ 3^2 

Hence, ( 1 ~ 1 > i_i 3 > 1- This is a contradiction since Ai > 0. Hence, 

there exists I > 1 such that p/^^(l) < 0. 
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Observe that (Ai (l — ■^)) = 1 for all I > 0. Hence, with same line 
of arguments as above, it can be shown that there exists I > 1 such that 

P/‘' Al (1 - si)) < 0- 

Now from (C.4) and (C.6) it is easily seen that > 0 for all I > 0. 

From the same relations we also observe that (Ai (l — ■^)) = 1 > 0 
for all I > 0. Combining the two we have 


(C.IO) ^max ^0, Ai > 0 

Further, observe that P{ '(Ai) = 0. Hence, there must exist at least one 

/'o\ 

root of Pj (a) in the following range 


(C.ll) 


cx G 


^max 




,Ai 


(2,) (2^ 

Let denote the minimum root of P| (a) in the above range. Therefore, 
in the range 


(C.12) 


cx G 


^max 




, min 



f 2') 

we must have P| (a) > 0. (Note that the right limit can be combined 

i2') (2) 

with 1 because of the minimality of ). Putting I = a = r\ in (C.6) 
we observe that p!^'^ < 0- Hence, using the same line arguments we 

conclude that in the range 


(C.13) a G ^max ^0, Ai ^ > ™fo (l; ^2^^) 


both P^\a) 


,{ 2 ) 


( 2 '\ (o') 

(a) > 0, where denotes the minimum root of P 2 (a) 


( 2 ) 

in the range defined in (C.12). Therefore by (C.6) we also have P( {a) > 
Pf (a) > 0 in the above range. Repeating the same argument again for 


Pg^^ (a) we find that PY'’ (a) > P 2 ^'' (a) > Pg^'' (a) > 0 holds in the range 


5 ( 2 ), 


>(2), 


>(2), 


(C.14) 


a G ^max 


0,Ai 


^ A 2 


, min ^1, 
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where denotes the minimum root of ' (a) in the range defined in (C.13). 

Trivially, we have (a) > 0 in the range defined in (C.14). Now from (C.5) 
we have P^^^O) = — AiA^^ < 0. Also, from definition of r® we know that 
P^‘^\r^'^) = 0. Now, by putting a = and / = 1 in (C.6) we obtain 



f2) 

Again, by putting I = 2 and a = rg in (C.7) and using the above we obtain 
Therefore, there exists at least one root of P 2 ^\a) in the 

interval ( ]. Denote the maximum of all such roots to be Hence, 

in the range 


(C.15) 


OL G 



, min 



we have P^^^ (a) > P^^^ (a) > 0 along with P^^^ (a) > P^^^ (a) > P^^^ (a) > 0. 
Again from (C.5) we observe that P 3 ^^(r f'x 0. Further, putting I = 3 and 
a = in (C.7) we obtain P 3 ^^(r 3 ^^) > 0. Thus, there must be at least one 
root of P 3 ^^(a) in the range ( ]. Let rg^^ denote the maximum root 

in the interval. Hence, in the interval 


(C.16) 


Of G 



, min 



we have (a) > P^^^ (a) > Pg^^ (a) > 0 along with pj^^^ (a) > P^^^ (a) > 
pj^\a) > 0. Similarly, from (C.5) we have P 4 ^^(r 3 ^^) < 0 and from (C.6) we 
have P^^ (rf^) > 0. Thus, there must be at least one root of P 4 ^^(a) in the 
range ( ]. Denote the maximum of all such roots by Hence, in 

the interval 


(C.17) 


Of G 



, min 



we have P 4 ^^(a) > P 2 ^^(a) > P 3 ^^(a) > P 4 ^^(a) > 0 and Pi‘^\a) > 

P^^\a) > P^^\a) > 0 . 
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Using the same line of arguments as above the following inductive hy¬ 
pothesis can be proved: If, for k > 0, P^^\a) > P 2 ^\a)... > > 0 

and Pl^\a) > P 2 ‘^\a)... > > 0 hold in the range 


(C.18) 


ot G 


max 1 Ai ( 1 - ^ j j , min (l,r 


,( 2 ) 

3+3fc 


then > P^^^a )... > (a) > 0 and P^'^^a) > p!^\a) 

-P 3 + 3 (fc+i) (®) A 0 hold in the range 


> 


(C.19) 


cx G 


max r 


zii 

4-h3(fc+l)’^l ^ 


,Ai 1- 


,min(l,rg3(;-^3) 


and the interval in (C.19) is included in the interval in (C.18). 
The decreasing sequence of compact intervals 


(C.20) 


max 


(rW 

1 '4-I-3A:’ 


Ai 



, min 



( 2 ) \ 
3+3fe J 


for k > 0 


eventually become strict subsets of the interval [0,1] as discussed in the 
beginning. Further, the intersection of all such compact intervals must be 
non-empty due to the Cantor’s intersection theorem. Hence, we have shown 
that there exists a G (0,1) such that the sequences |P;*'^^(q;), / G Z+| and 

^Pi^\a),l G Z+l are both positive non-increasing sequences of real num¬ 
bers in [0,1]. 

We now proceed to show that the above sequences satisfy property (P.3). 
Let lim/_).oo (a) = > 0 and lim/_).oo P/^\a) = ^2 > 0, where a is 

chosen such that both sequences become positive and non-increasing. Now, 
taking limit of (C.7) as / —)• oo we have 


(C.21) 


2 

E 

i=i 


11 

A, 



Now using the stability criterion and the fact that 0 < ^ 1,^2 E 1 we have 


1 1 

-r-1— 7 — > 1 

Al A2 


6 , ?2 


+ 


Ac 


>6>e 


d2 

2 
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with equality holding if and only if ^2 = 0. Further, we have 


Ai A2 - Ai A2 - 




d2 

2 


Hence, by multiplying both sides with we have 


6 , 66 

Ai A2 


> 66"^ 


— SI ?2 ! 


with equality if and only if 6=6 = 0- Again, since 6 ^ 1 '''^6 have 


^ 66 

Ai A2 - Ai A2 


> 66^^ > 


Hence, we have shown 


(C.22) 


iL + i6> 

Ai A 2 


SI S2 


with equality holding if and only if 6 = 6 = 0- Hence, for (C.21) to hold we 
must have 6=6 = 0- This proves (P.3) and thus completes the proof. □ 


APPENDIX D 

To prove Theorem 4.3, we first state the following lemma. We will write 
g < g' to mean that holds for all n € Z_|_ and j € J. 

Lemma D.l. If g < g' holds, for g, g' € , then u(t, g) < u(t,g') 

holds for all t >0. 

Proof. The proof is essentially the same as that of Lemma 3.3 of [9] and 
hence omitted. □ 

We define vlf\t,g) = Vn{t,g) = Y.jGjJjVn\t,g) 

for each n > 1 and j G H. Further, Vn\g) = '^k>n9k'^ ^n(g) = 

'^j£jljVn\g) for each n > 1 and j G J. 

Lemma D.2. //g G , then u(t,g) G for allt>0 and 


dvn{t,g) 


A 




i=i 


A, ) 


(D.l) 


dt 


for all n > 1. 
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In particular, 


(D.2) 


dt 


A 



A, ) 


Proof. Suppose that u(t, g) G holds for all t <t. Hence, ui(r, g) < 
oo and lini„_).oo ^ (t, g) = 0 for each j G J. Summing (4.4) first over all 
k > n and then over all j € H yields 


(D.3) 


dVn{t,g) 


dt 


= A 


t=T 



4-1 (eg 


M 

E 

i=i 


4 ') 


(h g) 


< oo, 


for all n > 1. Hence, for all sufficiently small h > 0, we have Vnir+h, g) < oo 
for all n > 1. This implies that u(r + /i, g) G for all sufficiently small 
h > 0. This fact along with g = u(0, g) G implies that u(t, g) G 
for all f > 0. Further, (D.l) can be obtained by summing (4.4) first over all 
k > n and then over all j € H □ 


Proof of Theorem 4.3. Clearly, Lemma D.l implies the following 


(D.4) u(t, min(g, P)) < u(t, g) < u(t, max(g, P)) 

Hence, to prove (4.16), it is sufficient to show that the convergence holds for 
g > P and for g < P. 

We first need to check that for each such g, the quantity ui(t,g) (and 
hence also Vnit, g) for n > 1) is bounded uniformly in t. If g <P, then by 
Lemma D.l we have u(t, g) < u(t, P) = P for all t > 0. Hence, ui(t, g) < 
Ul(P). 

On the other hand, if g > P, then by Lemma D.l u(t, g) > u(t, P) = P. 
Hence, we have 


(D.5) 


M 


E 


A, 


M 


> 


pU) 

E ' 


i=i 


A, 


= 1 


Thus, from (D.2) we have < 0. Hence, we have 0 < vi{t,g) < ui(g) 

for all t > 0. 
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Since the derivative of Un\t) is bounded for all j € J, the convergence 
u(t, g) —>■ P will follow from 


(D.6) ^ dt < oo, j E J,n > 1 

in the case g > P, and from 

(D.7) J - uW(t,g)^ dt <oo, j e J,n>l 

in the case g < P. Both the bounds can be shown similarly. We discuss the 
proof of (D.6). 

To prove (D.6) it is sufficient to show that 


(D.8) 



U^ri\t,g) - 


dt < oo, 


for all n > 1. We will use induction starting with n = 1. Using (D.2), we 
have 



A Jq dt 

= jiMs) - Mps))- 


Since the right hand side is bounded by a constant for all r, the integral on 
the left hand side must converge as r —>■ oo. 

Now assume that (D.6) holds for all n < P — 1. We have from (D.l) 
and (4.13) 
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^'L(0,g) - VL{T,g) 



dvL{t,g) 

dt 


dt 




U), 


uy\t,g) 


A, 


M ,\ 

n (4^ii(t,g)) j dt 




M 

n ( 4-1 (^>g)) 




dt 


M 


-a/ in(4-i(^’g 


0=1 


M 

n 

i=i 


P. 


U) 


L-l 


dt 


By the induction hypothesis, the last integral on the right hand side con¬ 
verges as r ^ oo. The left hand side also is uniformly bounded. Hence, the 
first integral on the left hand side also must converge as required. □ 
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